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quence of coefficients in the asymptotic expansion of the logarithmic deriva- 
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a measure, and also as the trace of the (n — l)th iterate of a linear operator. 
Applications of these sequences, and hence of the explicit solution provided, 
are found in quantum field theory as the number of Feynman diagrams of a 
certain type and order, in Brownian motion theory, and in combinatorics. 
Aequationes Mathematicae, 80, 291-318 (2010). 

Mathematics Subject Classification (2000). Primary 05A15; 
Secondary If B37,15A52,30Ef 5,33C15,35Qf 5,39B42,44A60,81T18. 



1. Introduction 

We study the sequence S'(ai,a2,a3) = (wn)J^Li defined by the self-convolutive 
recurrence 

n-l 

Un ^ (ain + a2)un-i + as^^UjUn-j, ui = I, (!) 

3 = i 

and show how to derive a closed-form solution as a Mellin transform, i.e. one of 
the form 

/■oo 

x"- V(a;) dx, (2) 

at least, in a reasonably broad special case (in general, there are extra terms, but we 
give an explicit expression for them). In other words, u„ is the (n— l)th moment of 
a measure. Clearly a representation such as ([2]) is useful if one wishes to compute u„ 
for large n — as it is nonrecursive — or study the asymptotics via Laplace's method 
[2] • Incidentally the representation additionally defines an analytic continuation of 
Un to Ren > 1. 
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Our study of such sequences arose from a particular example, namely the 
coefficients of the asymptotic expansion of the logarithmic derivative of the 
Airy function: 

Ai(z) ^ ^^A. (4^)3„/2 

with uo = — ^. As Ai'/Ai obeys a Riccati equation, the (it„) form a recurrence of 
the type given in ((Ij: 

n-l 

Un = (6n - 8)u„_i + ^ UjUn-y 

It was established in [10], thereby settling a long-standing question, that 
1 r a;"-4/3 

dx, n > 1, 



" 7r2 7o Ai(ia:2/3)2 +Bi(ia;2/3)2 

which is of the form © with ^(cc) = 7r-22:-i/3/ [Ai(ia:2/3)2 -|_ Bi(ix2/3)2] . This 
sequence and a related one described later are encountered in a number of ap- 
plications including the moments of the area function of a Brownian excursion 
[lOl [TT] . quantum field theory [5] (quite a fertile ground for these: more of this 
shortly) and combinatorics [3- For reference this sequence is 7^A062980 in the 
Online Encyclopedia of Integer Sequences (OEIS) [12] , and begins 

S'(6, -8, 1) (1, 5, 60, 1105, 27120, 828250, . . .)■ 

An obvious question to ask was whether this result can be generalised, and 
the answer is yes. A trivial example is S'(1,Q!,0) for which the nth term is T{n + 
a + l)/r(a -I- 2) and in Mellin form is 

foo a+l -X 

5(l,a,0)„^ / x-'^—-—-dx. 

Jo T{2 + a) 

However, it has a more interesting 'relative' (we shall say later in what sense) 



5'(1,-2,1)„ = / ,3 — jdx, 

where Ei denotes the exponential integral function [T]. It turns out that this se- 
quence has been studied extensively in combinatorics and algebra [7], because it 
is the number of connected, or indecomposable, permutations of [1 . . .n], i.e. those 
not fixing [1 . . . j] for < j < n. For reference this sequence is #A003319 in the 
OEIS and begins 

S{1, -2, 1) = (1, 1, 3, 13, 71, 461, 3447, . . .). 
The integral representation is, we think, a new result. 
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There are more known cases in this category, in fact, and here is another 
resuh from combinatorics: 

where 

oo r+- 

This sequence (#A000698 in the OEIS) begins 

3(2, -3, 1) = (1, 2, 10, 74, 706, 8162, 110410, . . .) 

and is the nmubcr of nonisomorphic connected Feynman diagrams of order 2(n+l) 
arising in a simphfied model of quantum electrodynamics (QED) [^5 . On multipli- 
cation term-by-term by 2n — 1 the following sequence (#A005416) is generated, 

(1,6,50,518,6354,89782,...) 

which is the number of 'vertex graphs' of order 2n arising in the QED perturbation 
series for the electron magnetic moment [4] . The above expression as the moment 
of a distribution has been stated by R. Groux [T9| though the result and derivation 
do not seem to have been formally published. This is one of the very few cases 
that have been solved. 

A sequence that bears intriguing similarity to that one is 

/■°° 2a;"~^e^/^ 
Jo G2{^x)^ +Tr^G3{^xy 

where 



G,{x)^2-G^{x)lnx + }_^^^:^, G,{x) . ^ , 

here Cr — (^^' ) is the rth Catalan number and 

dr = + 21n2 + ^ + 2E;.i (} - 

where 7 = 0.5772 ... is Euler's constant. This sequence (#A005412 in the OEIS) 
begins 

5'(2, -2, 1) = (1, 3, 18, 153, 1638, 20898, . . .) 

and is the number of Feynman diagrams with proper self-energies arising in QED 
[5] ; the integral representation is apparently new. 

Another example which generates a neater result is this, invoking the modified 
Bessel functions. The conciseness of its form is entirely due to the familiarity of 
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the functions involved, though, and in fact it is almost identical in construction to 
the previous one: 

2x"-2 dx 
Ko{\xr+7:^Ioi\xy 

2x"~2 dx 
K,Cixr+7r^h{\xr- 
The first few terms are 

S'(2,-4, 1) = (1,1,4,25,208,2146,26368,...) 
S'(2,0,-l) = (1,3,12,63,432,3798,41472,...). 

The first of these occurs in QED as the number of Feynman diagrams with exact 
propagators ^ and is #A005411 in the OEIS, though the expression as an integral 
is, again, apparently new. We are not sure if the second expression, which is not 
in the OEIS at the time of writing, has an application in the same field. 

This paper is organised as follows: the next section gives a complete exposition 
of the methods, and the one after that is devoted to particular examples. Notation 
is standard throughout; N denotes the set of natural numbers 1,2,.. .; — N denotes 
their negatives; Nq,— No denote the same sets including zero. We use F {x) to denote 
the digamma function, T' (x) /T{x); note that for m G Nq, we have F{1 + m) — 
-1 + E" 1 7 and F{\ + rn) = -7 - 2 ln2 + where 7 = 0.5772 ... is 

Euler's constant and the sums are disregarded for to = 0. 



5(2,-4,l)„ = / 



5(2,0,-l)„ = / 
Jo 



2. Methods 

2.1. Connection with Riccati equation, and solution as integral 

The use of the Riccati equation is fundamental, and we argue as follows. The 
sequence is generated from the asymptotic expansion of a function 0, i.e. 

-^W^E^^f^' Iarg^|<7r, (3) 

where we have added in the 0th term defined in the obvious way as un — — 7 — 
(if ai + a2 — 0, we consider the limit a\ -\- 0.1 — >■ 0). Note however that when 
subsequently we give formulas for u„, they will not necessarily be valid for n = 0. 
The function </) is assumed to be single-valued in the complex plane cut along the 
negative real axis, and in the situation at hand this condition will be satisfied in 
all cases bar one (the 'algebraic case'), which we deal with separately. 

The recursion in the (u„) corresponds to the following Riccati equation 
TZ{f3i, (32, 1^3, (34) defined below ([4]), in the sense that if we find a function sat- 
isfying (lU that has an asymptotic expansion ([3]), then the coefficients will satisfy 
^ and thus be the sought sequence: 

0'(z) = /?i0(z)2 + (/32 + (33/z)<j>{z) + (4) 
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with coefficients related by 

a "3 a ai+a2 + 2a3 1 "2 a 01+^2 + ^3 

Pi = , P2 = ^ ^ Pa = IH , P4 = 7 ^ (5) 

q;i Q;i(ai+Q;2j ai Q;i(ai+a2j 

(again, tlie case ai + a2 = can be dealt with by a limiting argument). Suppose we 
have found a solution to this equation with the appropriate regularity properties 
(in particular, lima;_>oo (j^i^) = ai+a^ )• Then we only have to extract the (u„) by 
the Fourier integral: 



27ri J^;; Z 

where the contour C runs anticlockwise in a large circle from —oo — ie to — cx) + ie 
(but it is not closed). 

It is probably fair to say that this last step has been something of a stumbling- 
block to explorers, and is the reason that the results in this paper have lain undis- 
covered for so long. On the face of it, the right-hand side of ([3]) is not a meaningful 
function, so apparently it cannot be valid to multiply by z"^^ and then perform 
the integral in such a cavalier fashion. However, the asymptotic series may be 
written as a finite sum of TV terms plus a remainder, Rf^(z) say, and crucially this 
remainder is bounded by cn\z\~^~^ in {\z\ > p, \ argz| < tt} for some constant cat 
independent of z. The integral can then be performed around a contour of finite 
radius > p and the contour expanded to infinite radius, causing the remainder 
term to vanish. Finally N is arbitrary and therefore can be made as big as we 
please. Incidentally we consider that much of the difficulty is caused by the '~' 
notation, which makes asymptotic analysis look like a branch of witchcraft not 
possessing the rigour of Taylor series. In fact, there is a symmetry between the 
two notions: whereas an asymptotic series has a remainder term that is bounded 
for a given N for all z inside some annulus, the Taylor series has a remainder term 
that is bounded for a given z for all n exceeding some N. For our purposes it 
is the first one that is useful. It is worth emphasizing that we are assuming the 
asymptotic expansion ^ to hold uniformly in the whole cut plane. This is not a 
problem in the case considered here, but for extensions or generalisations of the 
methods shown the region of validity needs to be borne in mind. 

This contour integral can now be evaluated as the sum of three parts, as 
follows. We collapse the contour on to the negative real axis: the resulting integral 
is the first part. Secondly, in wrapping the contour round the branch cut we 
may generate a contribution from the origin if the integrand is not integrable 
there. Finally, along the way we may have picked up contributions from other 
singularities of 4> anywhere in the complex plane except the negative real axis and 
the origin. 

It is convenient to classify cases by reference to their singularities. 

Definition 1. The sequence (or, equivalently, the associated Riccati equation) is 
said to be simply represented if it obeys both the following conditions: 

(i) (f> has no singularities other than simple poles in the cut plane C \ (— oo, 0]; 
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(ii) limz_j.o z(f){z) exists and is finite. 

In the context of this paper there are very few sequences that are not simply 
represented, and we can and do deal with them individually. We subdivide the 
simply represented sequences as follows: 

Definition 2. A simply represented sequence (or Riccati equation) is said to be 
regular if it obeys both the following stronger conditions: 

(i*) (p has no singularities in the cut plane C \ (—00, 0]; 
(h*) lim2^o z(j){z) = 0. 

It is said to be quasiregular if it obeys (ii) but not (ii*). Otherwise it is said to be 
irregular. 

Incidentally all the examples mentioned in the Introduction are regular. In 
the regular case we can collapse the contour around the branch cut and write 
z — —X to give 

1 f°° 

Un = Tr^ [<l)]_{x)x^-^dx (6) 

where [4>\_ {x) denotes the jump in 4){z) across the branch cut at z — —x, i.e. 
lime_j.o+ [0(— a; + ie) — 4>{—x — ie)] . So for a regular sequence, we have now obtained 
a solution in the form ([2|) with 

^^[x)^^^[<j^]_{x). (7) 

Note that 4> can be recovered from [0] _ immediately through the Hilbert transform: 

, . . 1 pI^LM, ^ 

(p[z) — —-—^ / dx, Rez > (J, 

27ri Jo x + z 

with the definition elsewhere being provided by analytic continuation, paying at- 
tention to the need to cut the plane along the negative real axis. 
In the quasiregular case we have an extra term: 

1 r „_i , . r -r\. n = l 



„ . , [(j)] (x)x"'-^dx , , ^ 

where r° is the residue at the origin. If desired this can be merged into /i(a;) as a 
delta- function, —r° ■ S{x). 

If there are simple poles in the cut plane, located at z = Q and of residue s^, 
further terms are added: 
1 f°° 

^n^TT t*^]- (^)^""' - • 0""' ' E (8) 

Writing the contribution from the origin as —r° ■ 0"~^, interpreting 0° — 1, makes 
it look like the contributions arising from other poles. 

Incidentally, simple poles on the negative real axis make their presence felt 
when the integral is collapsed on to the branch cut. A direct calculation of the 
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expression \mi^^Q^[(f){z + c+ie) — (j){z + c— ie)] shows it to be equal to — 27ri(5(x — c), 
with X = —z (z R^, but a quicker route is to use the Hilbert transform: 

6(x-c) , 1 
ax = 

Z + X z + c 

showing that the right interpretation of a pole in R~ in the generating function 
is a delta-function in the measure /it, of strength equal to minus the residue of the 
pole. 

If we were to relax the requirement (ii) in Definition [l] to allow a singularity 
worse than a simple pole, but still of power order, at the origin, then we would 
have to deal with a contribution from there of the form 

i.e. the integral around a small loop. By stipulating that z™'(j>{z) be bounded at the 
origin for sufficiently large to, we ensure that a singularity of this type can affect 
only finitely many of the We then say that the sequence is almost simply 

represented. 

We have shown in general how to solve for u„, but now need to give an 
explicit expression for 0, so we attend to that next. 

2.2. Solving the Riccati equation 

First, we exclude the case ai = 0, as then the Riccati equation is purely algebraic. 
This is better treated on its own, which we do in 21 

The next special case that requires attention is as — 0. In that case the 
Riccati equation is a linear differential equation and (j) is essentially an exponential 
integral function [1] §5]. We deal with this as the first of the special cases in the 
next section. It is well known that this function is regular in C \ (— oo, 0] and has 
a singularity of power order at worst at the origin. 

Another observation is the 'scaling law' 

S(tai,ta2,ta3)n = ^"^^(ai, a2, "s)™- 

By solving for (— )"~^Wn instead of Un, we can assume that ai > 0. 

Having dealt with these cases we are now free to assume (3i ^ 0, and so we 
can write (l){z) = —w'{z)/j3i'w{z) in the usual way to obtain the linear differential 
equation 

w" - [h + Pzlz)w' + PiliiW - 0. (9) 

The solution of this is essentially a confluent hypergeometric function (CHGF; [I] 
§13]), and the appropriate solution is 

w(z) oc e-("'=/'''"t/(a,6, kz) 

where 

fc=— , a = -/3i = ^, 6=-/?3--l-— ; (10) 
ai ai ai 
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notice that by assumption k, a, b are finite, k > and a^O. The effect of 'scaling' 
(q.v.) on these parameters is fc t^^k, with a,b fixed. Note that there is an 
apparent difficulty when 6 = 0, as then e~("''/'')^ is then undefined. However, 
the effect of the exponential is simply to add a constant on to w' /w, and in the 
asymptotic expansion this only affects Uq which is formally infinite, but not of 
interest. Unlike some of the other special cases (notably ai = 0, which requires 
special treatment), this can be dealt with by simply letting & — )• in the final 
expression. The function U is given (for 6 ^ Z) by 

U{a, b, z) = -Iii-^M(a, 6, z) + ^^^^z^-''M{a -b + 1,2 -b,z) (11) 

1 (a — + 1) i (a) 

in which M{a, b, z), also often denoted iFi(a; b; z), is the Kummer function, given 
by the Taylor series X^^q {ff" TT ' ^^'^ regular solution of the ODE 

zM" + {b- z)M' -aM = 0. 

Note that 

U{a, b, z) = z^-^lJ(a - & + 1, 2 - 6, z), 
which is one of 'Kummer's transformations'. A contour integral representation of 
U (a, 6, 0) for Re 2; > is 

1 



U{a,b,z) = -—I e-'H^'-'il + ty-''-' dt (Rea>0) 
1 (a) 



or 

r(0 



U{a,b,z) = ^ ^ ' / e-'H''-\l + t)'>-''-'dt (a^N) 

where in the second expression the plane is cut from to +00 and —1 to —00 and 
the contour runs from +00 just below the cut, leftwards, round the origin and then 
back to +00 on the upper side of the cut. The asymptotic expansion of U is 

U{a,b,z)^z 2^ , |argz|<f7r, (12) 

n=0 

SO the behaviour of U is known around the contour C mentioned earlier. 

We now need to evaluate [(/)]_ (x). From (f){z) = —w'{z)/j3iw{z) we deduce 

, ._k j^)U{aAkz) k 

a U{a,b,kz) b' 

Let U®{a, 6, z) and U^{a, b, z) denote the values of the U (a, b, z) for arg(2;) = +7r 
and — TT respectively. Then for z e M~, 



, kW{U®,UQ} 



^-(.)-<A-(.) = -^ 

with W {, } denoting the Wronskian. Now 



kz 



^® («' ^' ^) = ^^r^^-TXTT^("' ^' ^) + ^v^i^'~TM{a - 6 + 1, 2 - 6, z) 
1 (a — + 1) ^yO') 
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and similarly for the lower branch U^. So the Wronskian of and can be 
dealt with easily enough using the identity 

W {M(a, b, z), z^-''M{a - b + 1,2 - b, z)} = {1 - 6)^"^, 

and the result emerges as, on writing x for —z, 



W{U®{a,b,z),U^{a,b,z)} 



-2'Kix 



r(a)r(a- 6 + 1)' 



Hence 



/x(x) 



1 



k{kx) 



1 



(13) 



2^. ,^j_v-. r(a+ l)r(a-5+ 1) U®{a,b,-kx)UQ[a,b,-kx)' 

A little more work is needed on the last part. This involves the analytic continu- 
ation of U to the left half-plane, which is: 



U{a,b,xe^'''') = e 



r(i-6) 



T{a-b+l) 

T^i6 i-b r(b-l) 
r(a) 



M{b-a,b,x) 

M(l-a,2 



b, x) 



(14) 



We can therefore write 



with 



Ur{x) = e 



U{a, b, xe^''') = Ur{x) ± i[//(a;) 

r(i-&) 



r(a- 6 + 1 



-M{b-a,b,x) 



- (cos7r6) ^''^^ JK ^-''M{1 -a,2-b,x) 



V{a) 



Ui{x) 



(sin7r6) 



r(&~i) 
r(a) 



M(l - a, 2 - 6, x) 



(15) 



(16) 



and the U®U® term from above is expressed as 

U®{a,b,-kx)U^{a,b,-kx) = Unikxf + Ui{kxf . 
Now that wc have solved for wc can be rigoroiis about its singularities. 

Lemma 1. Suppose a\ ^ 0. The only singularities of (j) in the cut plane, if any, 
are simple poles of residue ai/as; unless as = 0, in which case there are none. 

Proof. Other than at the origin, the CHGF can have no singularities and nor can 
it have zeros of order > 1, as either would violate the differential equation that 
it satisfies. Hence the only singularities that w' /w can have, except at the origin, 
are simple poles of residue 1. So the residues of (p are always unless /3i = 

in which case (f) is an exponential integral and is free from singularities except 
possibly at the origin. □ 

Now let us be more precise about the singularity of <f) at the origin. By 
examining U (a, b, z) there we can deduce: 
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Lemma 2. Notation as above. Writing r° = lim^-^o we have: 

• Case a ^ —No, a — b + I ^ -Nq. If b > 1 then r° = (1 — b)/a; otherwise, 
r° = 0. 

• Case a e -No,, a - 6 + 1 ^ -Nq. Then r° = 0. 

• Case a ^ -No,, a-b+le -No. Then r° = (1 - b)/a. 

• Case a G -No, a-b+1 e -No (so b e Z). If b < 1 then r° = (1 - 6)/a; 
otherwise r° — 0. 

Proof. The first case can be obtained from (jlip , examining which of the two halves 
predominates when z — >• 0; see also [U §13]. For a G —No, U{a, b, z) is a polynomial. 
The second and fourth cases follow from direct examination of the coefficients. The 
third follows from Kummer's transformation to exchange a for a — 6 + 1, picking 
up a factor of z^^^ in the process. □ 

Incidentally if we assume that 4'{z) looks like c/ z at the origin then we can 
deduce immediately that c = or c = (1 — 5)/a directly from the Riccati equation; 
but choosing the right one is not quite so obvious. 

We have arrived at: 

Theorem 1. If ai > and ^ Q then £'(01,02,0:3) is simply represented, and 
hence solved fey Q. D 

Somewhat paradoxically, the ostensibly trivial case oi = has pecularities of 
its own that do not always permit a simple representation. We leave this case until 
the end of the paper because it bears resemblances to the general case that will 
only become apparent when a good selection of special cases have been exposed. 

Now we have to compute all the various bits. The following result summarises 
what we have found about ^i{x) and shows it to be a measure, up to a prefactor 
(remember that A; > by assumption): 

Theorem 2. The function /i in ^ is given by \13[ 15\ 16\i . After extraction of a 
constant factor it is positive for all x > 0. □ 

Now we turn to the question of zeros of U in the cut plane, which give rise to 
simple poles in (p. Denote by A the change in argument, divided by 2n, of U{a, 6, z) 
as z passes along each of the following paths: Ac, around the loop C anticlockwise 
from —00 to —00 as previously defined; A_|_, along the path from —00 to above 
the branch cut; Ao, for the clockwise loop around the origin; A_, along the path 
from to —CXI below the cut. Then by the Argument Principle, the number of 
zeros of U in the cut plane, counted according to multiplicity, is 

#[/ = Ac + A+ + Ao + A_ 

(remember that U has no poles in there). As the properties of U{a, b, z) on C and 
near the origin are well documented [TJ §13], we can track arg U very easily provided 
we can find A+ and A_, which we expect to be able to do from examination of 
Uii{x) and Ui{x), in certain cases at least. Note that A+ = A_ by symmetry. 

Lemma 3. The number of zeros of U{a, b, z) in the cut plane is as follows. 
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(i) // -2m - 1 < a < -2m + 1, m G No, and h < a+l, then #[/ = 2m. 

(ii) Ifa = b>l then #[/ = 0. 



Proof. In all cases one has Ac = —a, by (112]) . 

(i) We deal with the case 6 < 1 first. Write a = —2m — 1 + ^, with < 5 < 2, 
thereby temporarily excluding the case when a is a negative odd integer. The 
assumptions imply positivity of M(l — a,2 — b,x) for a; > 0, so 

^^(^) = r(2-5)r(a) ^ -M(l-a,2-M)| > q, 1<5<2 

Note that C/(0) = r(l - 6)/r(a - 6 + 1) > 0. From ([H]), and having just deduced 
that Im U{a, b, z) does not change sign as z passes from ooe+^' to above the cut, 
we must have A+ = ^{6 ~ 1). So #U ^ -a + ^{S - 1) + + ^{S - 1) = 2m. The 
case when a is a negative odd integer is more subtle because the number of zeros 
changes by two on account of a pair of zeros approaching the negative real axis. 
As the complement of the cut plane (in the Riemann surface) is a closed set, it 
contains its limit points, and so the extra pair of zeros remains in that set as a 
approaches and touches a negative odd integer from below. Hence, for example, 
with a ~ —3 we have #J7 — 2 and not 4. (When a is a negative even integer, i.e. 
6—1, U{z) simply travels along the positive real axis without going through the 
origin, so there is no difficulty.) 

The remaining case to deal with is0<a<l,l<6<a + l. Then Ui is 
negative above the branch-cut, by the same reasoning as above. The behaviour at 
the origin is different, as now U{a, b, z) ^ '^T(a)^ z^^^; this implies Aq = 6—1, and 
it is easily seen that A+ = i(a - + 1). Thus ifU^O again, just as we found 
above for < a < 1 and & < 1. 

(ii) From the behaviour of U at 0, we have Aq = a — 1. To attend to A+ and A_ 
we consider the behaviour of U{a, a, z)e^'^^^^^°'\ first on the upper side of the cut. 
This goes from —0 to +00 as z runs from ooe"'"'^' to Oe"'"'^', and we need to ascertain 
which side of the origin it passes on the way. By the analytic continuation formula 
for U we have, again writing z — —x, x G M+, 

U{a, a, a;e+^>-"'(i-'') = -e-^r(l - a)e'^"'Af (0, 1, x) - (real), 

so its imaginary part is — 7re^^/r(a) < 0. Accordingly, A+ — +i. So — 
-a + i + (a - 1) + i = 0. □ 

2.3. Symmetries 

Two symmetries are worth mentioning. First, by Kummer's transformation we 
have 

lU'{a,b,z) a-b+1 ( 1 U' {a - b + 1,2 - b, z)\ 1-b 



aU{a,b,z) a \a - b + 1 U{a - b + 1,2 - b, z) 
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Accordingly 



5(ai,a2,a3)„ = ^^^^±^±2^ ^(ai, -4ai - az, 2ai + ^2 + ag). 

201+02 „ ^ I 

03 ' 

0, n>l 



(17) 



which we call the Kummer symmetry. Two sequences can thereby be solved 'for 
the price of one'; as an example, 

5(6, -8, 1)„ ^ 5 5(6, -16, 5)„ - 4 • 

Secondly, the set of Riccati equations TZ{. . .) as previously defined is closed 
under reciprocation, i.e. the transformation 4>{z) i~> 6/<j)[z), which effects the fol- 
lowing transformation, 

7^(/3l, /32, /33, /34) ^ 7^(-/34/0, -/32, -/33, -A^); 

here 9 ~ — Uq — — (ai + 0,2)^"^ to make sure that wi = 1 in the transformed 
sequence. This induces the following transformation on the space of sequences: 

5(ai, a2, 03) ^ 5(ai, — 2ai — a2, ai + "2 + aa)- 

Writing ^ for the reciprocal, we have for example 

5(1, a- 1,0)^ = 5(l,-l-a,a) 

5(6,-8,l)t = 5(6,-4,-1). 

The effect on the CHGF parameters is 

k'^ — k; a) — a — h; — — &, 

and so 2a^ — U^ = 2a — b; also sequence is self-reciprocal iff 6 = 0. Clearly (5^)^ = 5. 
Note carefully that the reciprocal of a regular sequence need not be regular. 

2.4. Asymptotics 

The expression as the nth moment of a measure makes it particularly easy to 
examine asymptotics. For example, if 

IJ,{x) ~ Ax^e'^"" {x 00) 

then 

POO 

/ x''-^^l{x)dx AT{n + iy)\-''-''. 

Jo 

From (|12I13P we have 

, , kCkx^'^-^e-^'-' 

Uix) ^ ; ; -, X — >■ 00, 

' r(a + l)r(a-5+ 1)' 
provided the Gamma terms do not vanish, allowing us to write down: 

Theorem 3. // 5 is simply represented, a ^ — N, a — b ^ —N, fc 7^ 0, then 

r(n + 2a-b) 
u„ ^ — — r , n 00 . 

r(a-h l)r(a-6-h l)A:"-i 

Further, the reciprocal sequence to S has the same behaviour at leading order. 
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Proof. Clear from the above. Any poles in (j) generate power order behaviour, which 
is asymptotically dominated by the Gamma function. □ 

When a e — N the function U{a, b, z) is a polynomial of degree —a and then 
Un rises only exponentially, rather than factorially, in n. When a — b & —N we 
obtain the same deduction from the Kummer symmetry. 



3. Special cases 

In this section we deal with a variety of cases that arise in the literature, and we 
also deal with all the examples that we presented in the Introduction. The list is 
not exhaustive, of course, and is simply designed to give the reader a feel for how 
the method works. 

3.1. Factorial class S'(l, a - 1, 0) 
The solution is clearly 

_ T{n + a) 
"""" r(l + a) 

which is understood to terminate whenever —a G N (and (/> is a polynomial of 
degree —a then). The recurrence is linear and the Riccati equation is simply 



^'{z) = (l + -) 0(^) - - 



(where we temporarily restrict attention to a > 0) . The Riccati equation is solved 

by 

(j)(z) = —ze^'Eaiz) 
a 

with Ea denoting the exponential integral function (TJ §5]. Using the analytic 
continuation of to the left half-plane, we can write down 

and thereby ascertain that the (u„) are moments of the Gamma distribution with 
shape parameter a (note now that the result becomes valid for a > —1). 

When a < — 1 we obtain a sequence that is not simply represented, but almost 
simply represented: 4>{z) is of order at the origin and z(j){z) is not bounded. 
In the special case where a — —m, m S N, the sequence 'terminates' at and 
subsequent terms are all zero, and then (j){z) is of the form z^™ multiplied by a 
polynomial. 

This seems like a very heavy-handed way to prove an obvious result (and 
even more so if the CHGF is invoked) but, as we are about to see, we have not yet 
finished with the exponential integral function. 
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3.2. Recipro-factorial class S'(l, 1 — a, a) 

By reciprocating the previous sequence we can obtain more results. Indeed, 

Hz) 



azEa{z) 

or if we take the CHGF route, we have TZ{—a, —1, —a, 0) and a = b = a, k = 1, 
which gives 

wiz) = e-''U{a,a,z) = z^-°'Ea{z). 

Then 

Coo ^ \ 

r(l-a) + (cos7ra)xi-"y ^ 
^r + l-a r! J 

CO ^ J. 

UAx) = -(sin7ra)a;i-"e-^ V ^. 

r-O 

By Lemma [21 there are no zeros of Ea in the cut plane for a > — 1. (The case 
a = is unimportant; it is apparently singular because wo = ^, but this can be 
removed by subtracting the constant — ^ from (f>{z) before allowing a — > 0, and 
then one obtains 0(z) = —e^Ei{z), with fi{x) — e^^.) 

One subcase that simplifies a little is when a = ^ + m, m €E Nq. Then, the 
cos na term drops out and we are left with 

Ur{x) = e--r(i - m), Uj{x) = {-r+^x^-"^e-^ V ^. 

r— ^ 

For example, m — gives S'(l, — |, i) which on scaling by 2 gives the known result 
for 5(2, —3, 1) stated in the Introduction. (There is no singularity at the origin for 
TO = 0. For m > 1 one has r° = ^7^'" .) 

— l+2m 

When a € N the expressions need reworking as the r(l — a) term and the 
summation are both singular. The resulting expressions are 

Ur{x) ^ / ^ e-"Ei^(a:), Ui{x) ^ -^^^e-\ a = m e N, 
(to — Ij! (to — Ij! 

with 

Eim(a;) = Ina; — F{m) + (to — 1)! 



(r + 1 — m]r\ 

Now (/) is singular at the origin and, as Em{0) = rri > 1, a simple pole of 

residue is generated. In the special case m — I there is no extra contribution 
and we just have 



Ur{x) =e-^-l-\nx-J2^) = -e""Ei(x), C//(x) 



r=l 
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(note Eii = Ei), which gives 



5(1,-2,1): fi{x) = 



Ei(.T)2 + 7r2 ' 

as stated earher. More generally for a = m G N we have 



These sequences are classified in the OEIS (though the closed-form solutions are 
not; they are purely expressed through the asymptotic expansion of the Riccati 
equation) but except for m = 1 they do not seem to have any special combinatorial 
significance. 

Finally for a E —No we have another important case, as w{z) is then a 
polynomial x e~^. We deal with this later when considering a larger class that we 
are calling the Laguerre class. 

3.3. Airy cases S'(6, -8, 1) and 5(6, -4, -1) 

For 5(6, —8, 1) we have the Riccati equation as 7^(— i,0, — -i, ^), so the relevant 
CHGF solution is parametrised bya = i, h ^ i, fc = i. The solution is thereby 
given via 

w{z) = e-^/i2c/(i, i, iz) = 22/331/6^1/2 Ai(iz2/3) 

or equivalently 

Z-l/3Ai-(lz2/3) 

^^'^ ^ Ai(i.^/3) 

Now for X £ ]R+ we have 

i7ffi(i,i,_i2;) = e-^/i222/33i/V/2Ai(ic^x2/3) 
t/e(i^i^_i^) = 6-^/1222/331/ V/2Ai(icJx2/3), 

with (jj — e2'^'/3 and ZJ = e-2'^i/3 Using the identities 

Ai(e±2W32;) = ie±'^'/3(Ai(x) TiBi(a;)) 

we find 

;7©( 1 , 1 , _ i^)f/e( 1 _ 1 ^ _ 1^) ^ 2-2/331/3^6-^/6 [Ai(ia;2/3)2 + Bi(ia;2/3)2] 
whence by (|13l) we obtain 

^(6=-M): M-) = Ai(i.2;;^;Bi(1.2/3). - (19) 



as claimed earlier. The sequence is regular because Ai{jz^^^) has no zeros for 
2 ' 



I argz| < Itt (a known result il ). 
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The sequence S'(6,— 4, — 1) is its reciprocal, and the corresponding Riccati 
equation is 7?.(i, 0, ^, — ^) whose sohition is 



(^) 



zl/3Ai(lz2/3) 



4Ai'(iz2/3) 
which is generated from 

w{z) = e-^/i2[/(_ 1 , -1, iz) = -7ri/224/33-i/6Ai'(iz2/3), 

Hence by (IT3l) we obtain, as Ai' also has no zeros in the cut plane, 

^(6,-4,-1): Mx) = ^,(,J;^:^,^,,3). ^ (20) 

which begins 

S'(6, -4, -1) = (1, 7, 84, 1463, 33936, 990542, . . .)• 

The application of this is that it determines the moments of the area of the absolute 
value of a Brownian bridge [HI [TT] . 

As a brief demonstration of the asymptotic analysis that we derived earlier, 
we have the following table. For both sequences we have u„ ^ 6"(n — l)!/27r. 



n 


ln5(6,-8,l)„ 


ln5(6,-4,-l)„ 


Asymptotic 


1 


0.000 


0.000 


-0.046 


2 


1.609 


1.946 


1.746 


3 


4.094 


4.431 


4.231 


4 


7.008 


7.288 


7.121 


5 


10.208 


10.432 


10.299 


6 


13.627 


13.806 


13.700 


7 


17.224 


17.369 


17.284 


8 


20.972 


21.092 


21.021 


9 


24.850 


24.952 


24.893 


10 


28.845 


28.933 


28.882 



In this case the relative error in estimating u„ is ca. 5% for n > 10. 
3.4. Bessel class S'(l, 2v - 2, ^ - ly), < v < ^ 

This is a generalisation of the Airy example that we have just treated; they have 
in common the property that (32 = 0, or equivalently = —^{ai +0:2). The 
Riccati equation is TZ{f3, 0, 2/3, —1/4/3) with (3 = 1'—^, and the CHGF parameters 
are a = ^ — ly, b = 2a. The Airy cases taken above are given by v = ^, |. The 
sequence is irregular for ly outside the range [— |, |], and by Lemma |3] there are 
no zeros in the cut plane for < 1/ < |. The class is closed under the Kummer 
transformation, which simply effects the transformation 1— >■ — ^, so we need only 
analyse < < |, which is regular. The Bessel class is also self- reciprocal: the 
reciprocal sequence is parametrised hy ly^ = 1 — Note that the reciprocal may not 
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preserve regularity, e.g. u = —1 gives a quasiregular sequence but i/ = I — (—1) = 2 
gives an irregular one. 

In the Bessel case, the solution is 

w{z) = e-''^U{\ 2v, z) - t:-^'^z''K,{\z) 

(note that K^, = K^^), and 

1 K,-i{\z) 

The self-reciprocity can be seen from this solution, as (z^K^)' = —z'^K^^i and 
(z-'K,)' = -z-^K,+,. 

The cases ly — 0,1 may therefore be handled together and can be treated 
easily enough using the expansions of K^, Ki and their properties. We have 

K,{\z) Koi^z) 
4>\z) — 71— or -f— 

respectively (these can also be obtained from the relations K'^ — —Ki, (zKi)' = 
—zKq). The analytic continuation of Kq{z) round to the left half-plane is 

K®{-x) = Ko{x) ~ TTiloix), K^i-x) = Ko{x) + 7ri/o(a;), 

so we obtain 

On scaling by 2 one obtains the results mentioned in the Introduction. 

When v is not an integer a direct appeal to the CHGF has to be made, and 
we evaluate Un and Uj using the result 

M(i + 1^,1 + x) = rr(l + v)x-''e''''^I^{\x). 

The duplication and reflection formulae for the Gamma function allow all the T 
terms to be cleared up, and the result is 

Ur{x) = TT-i/^^-'e-^/^ (iC^(ia;) +7r(sin7ri.)/^(ix)) 
Ui{x) = -7ri/2(cos7ri/)a;''e-^/2j^(i^)^ 

where we have also used the expression for in terms of I±^, 

K,{z) = -^^{l-.{x) - I,{x)). 

2 sm t:v ' 

The above expression for Ur, Ui is also valid for integer v. 

By (|13p we obtain 

(i — U)^^{C0STTI')/X 

2- -I- t-y^J - K^i^xy + TT^I.i^x)^ + 2Trismniy) K,i^x)I,i^x)' 

(22) 



S{l,2iy-2,^-iy) : n{x) 
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As expected, the Airy case can be obtained from this by setting v — \^\■ Now 
the cases v = \ cause minor difficuhies in ([22]) . However, \i{x) just boils down 
to in the case = i, and so is not of particular interest given that we have 
already dealt with the sequence S'(l,— 1,0). When = | we have S'(l,l,— 1) 
which is simply the powers of 2. This is because 0(z) has a pole on the negative 
real axis (at z = —1), and indeed the denominator of ^2f2\ completes as a square 
that vanishes at a; = 2, so is to be interpreted as 8{x — 2). 

3.5. Miscellaneous 

The sequence S{2, —2, 1), which is self-reciprocal, has a — ^, b = 0, k — ^. From 
Lemmas [2] and [3l the problem is regular. The expression for Uji{x) is singular at 
6 = by virtue of the Gamma functions, but the singularity is removable and 
application of L'Hopital's rule yields 



Ur{x) 

with 
and 



2 lnx^ r(^+r)x'-+i ly-£(i_LLZ^ 
~^r(2 + r) r! ^ ^ ^ r(2 + r) ' 

dr = r(2 + r) - r(i + r) + r(l + r); 



r] , x'^^ 



r^r(2 + r) r\- 

After some minor manipulations one ends up with the result given in the Intro- 
duction. 

The related sequence (#A005413) 

(1,7,72,891,12672,202770,...), 

which is obtained by ti„ = (n — l)(u,i + 2rtu„_i), n > 2, is thereby solved too, and 
in QED is the number of graphs with 'proper vertices' [5]. 

3.6. Laguerre class 5(— 1, 2 + a, to), S'(— 1, 2 — a, a + m), m G N 

The second of these cases is the Kummer transformation of the first, so we only 
have to deal with the first. In this case k — —1, a = —to, b = \ -\- a and then the 
CHGF is simply a polynomial. (We have chosen to reverse the signs and have k 
negative so as to keep the branch cut along the negative real axis.) The solution 
is given by 

w{z) = e-'""/(i+")C/(-TO, 1 + a, -z) 
which is expressible in terms of the Laguerre polynomial defined by 

Lt\z)^^-^U{-m,l + a,z) = ^ i e-^*t-'"-i(l +<)'"+" dt. 
to! 27ri 7(0+) 

One therefore has 

Lt^'{~z) , 1 



{■A 



mLm\-z) 
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which admits a partial fraction expansion 

1 1^1 

j=i 'J 

Intriguingly, the sequence therefore obeys a linear recurrence relation of order m, 
with constant coefficients. At this point one might jump to the conclusion that as 
the zeros of the Laguerre polynomials are well-studied, everything must be plain 
sailing. That is not so. 

When a > — 1 the zeros of Lm \x) are known to be real and positive, as they 
are orthogonal polynomials of the Gamma distribution (density oc a;"e~^). The 
theory of these is extensive. For example, the zeros of Lm (x) interlace those of 
Lm]_i{x); and asymptotic results are known for large m. The sequence is therefore 
regular in this case, and one has, on writing for the zeros of Lm\x), 

5(-l,2 + a,m)„ = -Ver' («>0). (23) 

i—l 

So the sequence has a Mehin representation as a 'comb' of deha-functions 

K^) = — y''5(a;-C*)- 

m ^ — ' 

i=l 

This is a discrete probability measure bearing some resemblance to the Gaussian 
quadrature formula of the associated Gamma distribution: the locations of the 
delta-functions are the same, but the strengths of the delta-functions are all equal 
to — , unlike the ordinates of the Gaussian quadrature formula. As an example, 
a — 0, m — 2 gives 

2, 2)„ = i((2 + ^/2)"-l + (2 - V2)"-i), 

so the delta-functions are situated at x = 2 ± a/2- It is readily verified that u„ = 
4u„-i - 2u„_2. 

But when a < — 1 the connection with Gaussian quadratures disappears, and 
now the zeros need not be real. In general, these sequences are not regular. The 
relation (1231) is still valid but there is in general no representation as a measure. 
As examples, with a = — 1 and m — 3 we get 

S{-1, 1,3)„ - i((3 V3)"-i + (3 - V3)"-i + 0«-i), 

which is quasiregular. On the other hand a = —4, m = 3 gives 

S{-1, -2, 3) = (1,-1,-1, 3, 3, -21, 27, 27, -117, 27, 459, -837, . . .) 

which is irregular (and in appearance too) and solved by (j23p with roots given as 

- vV2- 1 - vV2 + 1 - 1 
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these being the roots of + + 6x + 6 = 0. The associated hnear recurrence is 
Un = —Zun-i — 6u„_2 — 6u„_3. Notc incidentaUy that up to scaling by —1 this is 
S{1, 2, —3) which is in the recipro-factorial class. 

One thing that does carry over from the theory of Gaussian quadratures 
irrespective of a is the connection with Jacobi, or tridiagonal, matrices [B]. By 
exhibiting an m-dimensional tridiagonal matrix T whose characteristic polynomial 
is the monic Laguerre polynomial, we can simply write (j23p as 

7.„ = -tr[T""i], (24) 

TO 

where tr denotes the trace and is the identity matrix. An explicit expression 
for T is 

" 1 + a 1 

1(1 + a) 3 + a 
y= 2(2 + a) 

m'{m' + a) 2m' + 1 + a 

where to' = m — 1. 



1 

5 + a 



(25) 



4. Algebraic class 5'(0, a, 1) 

We round off this paper by treating what is on the face of it a rather trivial 
special case: when ai — 0, the Riccati equation reduces to a quadratic (algrcbraic) 
equation. It turns out that this case is distinct, one reason being that different 
types of singularities are encountered from the general case that we solved using 
the CHGF. 

This case is definitely of practical interest as it crops up frequently in com- 
binatorics, most obviously through the well-known Catalan numbers, 

n 

Cn — 5(0, 0, l)n-|-i; C'n+l — CjCn-j- 

j=o 

We now proceed to solve this case for all a. The Riccati equation is solved as 

(1 + 2/a)z + a - Jz'^ + 2(2 + a)z + 

H^) = 

which needs careful interpretation, as follows. 

• a > 0. The plane is cut between z = ~2 ~ a ± 2^/1 + a and the square root 
in the above equation is taken as real and positive when z is. At z = the 
numerator of the above expression for (p is a — \fo? = 0, so ^ is regular there. 
By considering the jump in </) across the branch cut, we identify the series as 
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regular and with a measure that is bounded: 

( J-x^ + 2(2 + a)x - a2 

I 0, otherwise. 

• a = (Catalan). We can let a — >■ in the above case without vitiating the 
result, despite the singular nature of ^ (which is caused by uq = a^ai ~ 
and hence is of no consequence). This gives 



\/4 - X 

Ax) = \ l^-2|<2 (27) 

0, otherwise. 

— 1 < a < 0. The plane is cut as before with the same definition of the square 
root. But now a — \fa^ = 2q:, and so ^ has a simple pole at the origin of 
residue a. Thus a delta-function is created: 



{J-x^ + 2(2 + a)x - , ,, „ ^ 
' \x-{2 + a)\<2VT+^ (28) 
{—a)S{x), otherwise. 

a = —1. Here (j){z) = —1/z and so iJ.{x) = S{x). 

—2 < a < —1. The branch points are complex with negative real part; 
again the square root is taken to be real and positive whenever z is. Write 
a = + a\ and parametrise the complex variable on the branch cut by 
Cj/ = — 2 — a + 2\ay. We arrive at 



-— f {-CyT-^^l^dy, n>l. 

J-l 



The contribution from the pole at the origin takes care of the ui term. 

• a = —2. Here (j){z) = — (1 + y^^/4 + 1) / z. The even-numbered terms are the 
Catalan numbers, with alternating signs, and the odd terms vanish after the 
first: 

5(0, -2, 1)2„ = {-TCn-l = -- C (-V)»-ly/r^rf?/ 

J-l 

This can also be seen from a direct appeal to the Binomial expansion of 4'{z). 

• a < —2. This case is similar to— 2<a<— 1, but the branch points are in 

the right half-plane and the origin lies to the left of them, so <j) is now regular 
at the origin. We therefore obtain, as before, 



J-l 



■CyT-^Vl^dy, n>l. 



A convenient summary is: 



Theorem 4. The sequence 5(0, a, 1) is regular if a> 0, and quasiregular if —1 < 
a < 0. If a < —1, then it is not simply represented. □ 
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It appears that there is an interesting connection between some of these cases 
and the theory of random matrices, and we pursue this discussion in the 
last section of the paper. 

Trace formula. In aU cases for a > —1 the sequence is regular and admits a trace 
formula in the usual way. The other situations do not quite fit the pattern we have 
previously encountered. For a < — 1, the integral 

L 

bears a resemblance to a continuous version of the expression X^jlil^O)"^^ 
that we encountered in the irregular Laguerre case, and it is tempting to think of 
it as an integral over some 'pole density' However, the measure would have 

to be 

which is not real, let alone real and positive. So there is no obvious interpretation 
as a limiting distribution of poles arising from the Laguerre case. An open question 
is, therefore, how to interpret the complex density, possibly by trying to split it 
into real and imaginary parts and examining each separately. 



5. Conclusions and Discussion 

We have considered the sequence ll} and shown how it may be solved in closed 
form as the expression ([8|). In the case we have described as regular the solution 
reduces to 

We list some questions and applications that naturally arise from this work. 

Trace formula. It seems that the expression for m„ is the trace of the (n — l)th 
power of a linear operator. In the Laguerre case the spectrum of this operator 
is discrete and it admits a matrix representation. Can this representation as the 
trace of an iterated operator be extended to all the sequences? It seems reasonable 
to suppose that the operator has something to do with the physical system that 
gives rise to the sequence. 

Random matrices. It appears that there is an interesting connection between some 
of these cases (PM^5|) and the theory of random matrices. Let X he & {q x p) 
matrix whose entries are independent and Normally distributed with mean zero 
and variance 1. Consider now the empirical covariance matrix 

S = \ tXX', 
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in which ' denotes transpose, and let p, g — > oo with p/q held fixed. The limiting 
distribution of the eigenvalues of this {q x q) matrix is, provided p > q, 



^{\+-x){x-\_) 
, A± = (1 ± 



27r(g/p) X ' 

which up to a scaling is identical to p{x) in eq. (j26p . This result was derived rigor- 
ously by V. Marcenko and L. Pastur 12 and studied more recently by A. Sengupta 
and P. Mitra [TB]. 

When — 1 < a < 0, the connection still works. This is the case when p < q, 
so S is rank-deficient and must have at least q — p zero eigenvalues: hence the 
occurrence of a point probability mass at zero in (j28p . Of course the case p = q 
is covered by (l?7|) . which on changing variables x ^ x'^ states that the square 
roots of the eigenvalues are distributed according to the 'semicircle law' \/4 — x^/tt 
(0 < a; < 2), a famous result of E. Wigner [201I2T]. 

Both of these can be understood as limiting instances of the Laguerre class 
5(— l,am, m) as m ^ oo. The Laguerre polynomials also arise in the theory of 
random matrices, and clearly this is an avenue worth investigating. 



Numerical analysis: Finding zeros of U{a,b,z). When the sequence is irregular, 
(p has poles in the cut plane. So far the emphasis has been on finding a closed 
form solution to m„, which would necessitate counting them and finding their 
locations. However, if we use the methods in reverse, and compare the integral 
Jp x"~^^{x)dx with Un for a few values of n, we can find them. An example is 
adequate to demonstrate the general case, and we take the Bessel case with z/ = 2. 
We have a = — |, 6 = —3, and there is no contribution from the origin. By Lemma 
[3] there are two zeros to find. We therefore have 

where K2{^Cj) — 0- Substituting n — 1 simply confirms that we were right about 
the number of zeros, as the integral (without the prefactor) is equal to ^. For 
higher n we evaluate the integral numerically, which presents no difficulties as the 
integrand is smooth and unimodal. We arrive at the following table: 



n 




Integral 


Difference, 


1 


1 


-0.333333 


1.333333 


2 


2.5 


-0.916997 


3.416997 


3 


5 


-2.77313 


7.77313 


4 


5.625 


-9.25324 


14.87824 


5 


-15 


-34.4717 


19.4717 


6 


-154.6875 


-145.797 


-8.89074 
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The zeros obviously occur as a complex-conjugate pair, re^'^ say. Therefore the 
difference is equal to 



u„ = |r"-i cos((n - l)(7r - 6)). 

We now employ Prony's method (see e.g. [T3]) to extract r,9. This consists in 
observing that (?;„) obeys a linear recurrence, here of order 2: 

V„ + CiU„_i + C2W„_2 = 0, 

with the coefficients identifiable easily enough from (w„) by linear regression, af- 
ter which the (— Cj) £^re the roots of the quadratic equation -t- ciz + C2 = 
(equivalently, ci = 2rcos0 and C2 = r^). Applying this to the above, we find 
ci = —5.12549, C2 = 7.305500, from which we deduce that the zeros of K2{\z) in 
the cut plane are located at 

C - 2.70287e±2-8i8i8i. 

A few points can be made about the above analysis. First, we have identified the 
complex zeros of K2{\z) despite evaluating the function at real values only, which 
has a certain appeal to it and is particularly useful if one does not have at one's 
disposal an algorithm for its evaluation at complex arguments. The method is 
applicable to higher orders of modified Bessel function, which have more zeros. 
In practice, we suggest that the best way to find zeros is probably to use the 
technique to obtain reasonable accuracy and then use Newton's method to polish 
them one by one. On the zeros of K2[\z), using the above estimate as starting- 
point, we find the positions of the roots (using the computer package MATLAB) to 
be in agreement with the above, to 6 s.f. This method has an important advantage 
once the number of zeros becomes reasonably large. Without such assistance, one 
would have to search for the zeros of a special function without a clear idea of 
where to find them (Newton's method not being globally convergent in general). 
The associated polynomial in Prony's method is of degree #J7, so for #J7 < 4, 
it can be factorised by radicals; for > 4, search algorithms have to be used 
to factorise it, but this is a reasonably standard problem |16j . and certainly much 
easier than finding the zeros of a special function. 

Number-theoretic properties. So far, the analysis of the recurrence has been done 
using complex variable theory. However, there are certain questions that are likely 
to be resolved more easily by other methods. For example, for how many n does 
Un = 5 {a question that is only interesting for the irregular series in which terms 
vary in sign, as in S'(— f,— 2,3) earlier); or for now many n does u„ = 3 (mod 
7)? Questions of this form are much easier to tackle using p-adic analysis [3j [l7] 
than complex analysis: one such theorem that is arises from local field theory is the 
theorem of Mahler & Lech, which states that if a sequence obeys a linear recurrence 
with constant coefficients, then the set of n for which m„ = (or some other 
specified number) is a finite combination of sets each of which is either periodic or 
finite. Obviously this theorem does not apply to the recurrence discussed herein, 
except for the sequences of the Laguerre class, which do obey a linear recurrence. 
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However, the hypergeonietric functions used in the paper do havep-adic analogues, 
even though the convergence properties are generally very different in the p-adic 
valuation (essentially because \n\\p — > as n — oo), so the use of p-adic analysis 
may be fruitful. 

Other recurrences. The methods used in this paper could possibly be applied to 
other self-convohitive recurrences of a more general type. It seems that quantum 
field theory yields a rich source of integer sequences obeying similar self-convolutive 
recurrences. For example, L. Molinari |14] gives some sequences related to the enu- 
meration of Feynman diagrams that arise from asymptotic expansion of solutions 
of first-order nonlinear differential equations that have a cubic nonlinearity, 

(we have altered his notation so that the expansion is in ascending powers of 
—z~^), and hence require further manipulation before they can be reduced to 
a Riccati equation. A recent article by Y. Pavlyukh and W. Hiibner [TS] solves 
for the generating function but does not give a final reduction to an expression 
resembling Q; hence further work is required in this direction. There is also the 
possibility of using more general types of differential equation such as the general 
hypergeometric equation, and possibly invoking the Painleve transcendents as a 
more general class of nonlinear differential equations. 

That these sequences arise in other branches of mathematics should make 
this work quite widely applicable. Aside from this there is an elegance about the 
result that makes the subject quite attractive, and we hope that this, combined 
with the applications to the problems from which the sequences arise, will give 
other researchers the enjoyment that we have had from investigating this topic. 
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